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^. Abstract 






In this paper, we study the dynamics of a random walker diffusing on a disor- 
dered one-dimensional lattice with random trappings. The distribution of es- 
cape probabilities is computed exactly for any strength of the disorder. These 
J3 ' probabilities do not display any multifractal properties contrary to previous 

numerical claims. The explanation for this apparent multifractal behavior is 
'^ , given, and our conclusion are supported by numerical calculations. These ex- 

^ \ act results are exploited to compute the large time asymptotics of the survival 

O \ probability (or the density) which is found to decay as exp[— Ct^/^log^'^(t)]. 

An exact lower bound for the density is found to decay in a similar way. 
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I. INTRODUCTION 

The dynamics of the survival probabihty of particles diffusing in the presence of traps is 
a rich problem which has been widely discussed in the physical and mathematical literature 
within the last two decades ||l|-§[. 

The simplest system is that of diffusing particles in the presence of perfect static traps 
]l|-^. This problem (that we will call the Donsker-Varadhan problem) has been solved using 
very different technics. The main results is that the density does not decay exponentially 
(as a simple mean-field argument would predict) but as, 

n{t) ~ exp -CdC^+^td+2 , (1) 

where c is the trapping site density. The physical interpretation in d dimension is that the 
process is dominated by particles standing in very large trap-free regions of linear size L 
(these regions have a probability of order exp(— cL'^)). In such a region, the density decays 
as exp(— t/L^). A saddle-point argument then leads to the result of Eq. (|ip, with the relevant 
regions being of typical size L ~ (t/c)*+^, at time t. 

In another class of models |^, the traps are allowed to move. When these traps undergo 
free diffusion the density of particles decays as, 

n(t) ~ exp -CdCt^ , (2) 

for d < 2, and decays exponentially for d > 2. This result holds in the case of static or 
diffusing particles |^. 

It would be interesting to introduce the effects of hopping disorder on the trapping 
process. Even without trapping, quenched disorder in the particle hopping probabilities is 
known to have very important effects on the diffusion and first return properties ^-^- In 
the case of symmetric hopping probabilities (wj,i+i = Wj+i.j) 0,0, anomalous diffusion is 
observed, with an exponent depending on the properties of the disorder. In the generic 
non symmetric case (see for a more precise criterion), as in the Sinai model P] where a 
particle diffuses in a random (Brownian) potential, the diffusion is dramatically suppressed, 
the particle being effectively trapped in deeper and deeper valleys of the potential as time 
goes on. 

In the present article, we study the dynamics of particles diffusing in a symmetric or non 
symmetric disorder, in the presence of a random finite trapping probability at each site. 

II. MODEL AND KNOWN RESULTS 

Consider a particle moving on a one-dimensional lattice with random barriers (or hopping 
probabilities) and random trapping probabilities. More precisely, a particle at site i has a 
probability Wi,i+i < 1/2 (respectively Wj^j-i < 1/2) to hop on site i + 1 (respectively i — 1), 
and a probability (1 — 7)(1 — Wi,j_i — Wj^j+i) to disappear (7 < 1). With residual probability 
7(1— W7j^j_i — Wj^j+i), it just stays on site i. The hopping probabilities can be taken symmetric 
(wj^j+i = ti'i+ii) or non symmetric, and will be chosen according to the following typical 
probability distribution: 



p{w) = 2^-^(1 - f3)w-'^e{w)e {1/2 - w), (3) 

where /? < 1 measures the quenched disorder strength. 

The case 7 = 1 (no trapping) has been extensively studied ^^- In the symmetric case 
I^J^, one observes anomalous diffusion, {x^it)) ~ t^*^, with 1/ depending continuously on /3. 
The return probability Ps{t) = {Pi.i{t)), that is the probability of being at site i at time t 
having started at site i, decays as Ps{t) ~ t^'^"^'^, where ds is the spectral dimension [§|,0- 
The Sinai model |^ describes the generic non symmetric case and displays logarithmically 
slow diffusion and other peculiar properties. 

In the presence of trapping (0 < 7 < 1), the problem has been studied essentially by 
numerical means |]9|JT0[|. In addition to Ps{t), one can define the normalized return probability 
P{t) as, 

Note that in order to keep the notations simple, it is understood that (. . .) involves an 
average over the disorder and the considered site i. The proper way of defining (a;^(t)) is 
now, 

only taking into account surviving particles. With these new definitions, (a;^(t)) ~ t^*^, with 
2i/ Ri 1.25, seemingly independent of 7 and the disorder strength (3 |]10|- P{t) decays as a 



power-law, P{t) ~ t ", with a ~ 0.59, also independent of 7 and (3. Due to trapping, the 



survival return probability Psit) = {Pi,i{t)) decays much faster, and the authors of ||T0[ gave 
an heuristic argument leading to, 

logP,(t) ~ -Vi, (6) 

in qualitative agreement with numerical simulations. Moreover the probability distributions 
of quantities such as Pij{t) have been shown numerically to be very broad, leading to non 
self- averaging effects. 

Some of the peculiar properties of this model have been related to the possible existence 
of multifractal distributions for quantities such as the escape probability (see below) , in 



analogy |Tl| to what has been observed for the Sinai model [§[. 

For instance, let us consider the probability Gi^i+i(t) that a particle makes a first passage 
from a site i to a site i -\- 1 in t steps. This obeys the following master equation [Q : 

Gi,i+i(t) = wu+iSt,i + wu^iGi^ij,+i{t) + 7(1 - wu+i - Wi^i^i)Gi^i+i{t - 1), (7) 

with the boundary conditions, 

G'o,i(t) = woJt,i + 7(1 - wo^i)Go,i{t - 1), (8) 

and Gi^i+i{t) = 0, for i < 0. Now the escape probability is defined as. 



X, = X:G,,,+i(t), (9) 

t=o 

and the authors of P claim to have found a multifractal distribution for this quantity, after 
averaging over i and the quenched disorder. 

In the next section, we show analytically that this escape probability has in fact a well 
defined distribution and that the apparent multifractality can be ascribed to peculiar features 
of this distribution. In addition, we use some of these results in section IV to derive the 
exact large time behavior of Psif), by computing the properties of the survival probability 
distribution. We also give an exact bound on the density of particles, which fully confirms 
our result: 

\ogPs{t)^-t^'^\og^'\t). (10) 

Sections III.C and IV. B, where a rigorous method to analyse Lifshitz-like tails is in- 
troduced, are rather technical. The reader more interested in the physical consequences of 
these results could skip these technicalities and take the qualitative arguments given in the 
beginning of these sections for granted. 



III. ESCAPE PROBABILITY DISTRIBUTION: NON SYMMETRIC CASE 

A. Preliminaries 

Let us show that Xj+i can be simply evaluated from the knowledge of Xj, a calculation 
already appearing in |T^J^. We define the generating functions, 

+00 
G,,{z) = Y.z'G,,{t), (11) 

t=o 

with Xi = Gi^i^i{z = 1). The convolution theorem ensures that, 

Gi-i,i+i{z) = Gi^i^i{z)Gi^iJ^i{z). (12) 

Using the master equation Eq. (|^) this straightforwardly leads to the iterated map: 

Xi = 7- '- — ^ , i > 1, (13) 

1 -7(1 -pi -qi) -QiXi^i 



and. 



l-7(l-po) 



where Pi = Wj^j+i and g^ = Wi^i^i are random variable of identical distribution p given in 
Eq. (^. Note that these variables are independent only in the non symmetric model, as the 
extra constraint Qi = pj_i holds in the symmetric case. 

To the price of lengthier calculations, we have checked that the asymptotic behaviors 
of the observables of interest are not affected by the value of 7, provided trapping is not 



suppressed (7 < 1). This is in agreement with the numerical simulations performed in P,p!0 . 
From now on, we therefore restrict our analytic study to 7 = 0. 

In the case 7 = 0, the maximum possible value for x should satisfy, 

Pmax (-\r.\ 

J- O'lnax-^max 

with Prnax = ^max = 1/2, leading to Xmax = 1- 

B. Equation for the stationary distribution: non symmetric case 

As noticed above, the p^'s and g^'s are independent variable in the non symmetric case. 
Let us now assume that for large i the probability distribution of Xi exists and becomes 
independent of i, in contradiction with the numerical claim of 0. The stationarity of the 
distribution is exploited by expressing that the distribution of x, should be the same as that 
of Xi-i (at least for large i), leading to: 

H /-1/2 H/2 / p \ 

fix) = dy dp dqf{y)p{p)p{q)S x - . (16) 

h h h \ \ — qy ) 

This equation, though apparently very complicated, can be in fact exploited quite precisely. 
First consider x < 1/2, p = x{l — qy) is then always within the integration domain [0, 1/2], 
for any q G [0, 1/2] and any y G [0, 1]. We then find, 

f{x)=Io'dyJ^^'dqf{y)p{x{l~qy))p{q){l-qy), (17) 

= 21-/3(1 - /3)x-(^J,Uy!,'/'dqfiy) (1 - qy)''^piq), (18) 

that is a pure power-law behavior. 

For X > 1/2, imposing p = x(l — qy) G [0, 1/2] leads to new integration bounds: 

fix) = 22(1-/3) (1 _ p)^^-Pjl^^, rfy /(IV)-).- dqfiy)il - qy)''^q-^, (19) 

= 2^(1-/3) (1 - p)^x-^ Jl^_, dyftU-^ dqf{y)yP-\l - q)'-^q-^. (20) 

Introducing g{x) = x^f{x), one can differentiate twice this equation with respect to the 
variable x, leading to the following non local differential equation for g, valid for x G [1/2, 1]: 



g"ix) + 



3-/3^ (3 



X x'^ (2 — X 



-(i+,0) 



g'(x) = (1 - /?)2x/3-5 [2 - x-1] ''^" g{2- X-') . (21) 



The structure of this equation is quite unfamiliar as the LHS linear differential operator is 
determined by the value of ^f at a;' = 2 — x~^ in the RHS. 

g{x) = x^ f{x) being zero for x' < (that is x G [0, 1/2 [), this equation actually re- 
produces that g'{x) should be zero for x G [0, 1/2 [, in agreement with Eq. ([T8|) . Then, the 
knowledge of g{x') for x' G [0, l/2[ determinates the LHS on the interval x G [1/2, 2/3]. g is 
then determined by imposing that it is continuous at x = 1/2, and a consistency equation at 
X = 2/3 (see below). This procedure can be iterated. Consider uq = 0, and Un+i = 1/(2— «„) 



(that is Un = 2 — u~li), the knowledge of g on the interval [u„_i,m„[ and the continuity 
condition at x = u„ for g and its first derivative fully determine the function g on the next 
interval [un,Un+i[. Un can be exactly computed by induction leading to, 



l-Ur^ 



n + 1 



(22) 



As can be expected, u^. 



Xin 



1, when n -^ +oo. 



This recursion process and the form of Eq. (^) ensures that g is infinitely differentiable 

[. Moreover, if g is dn+i times differentiable at x = Un+i, it is at 



on the interval 



Pru Un+l 



least dn+i + 1 times differentiable at x = m„+2- This shows that the continuity condition 
for g and its first derivative suffices to determine g on the interval [un,Un+i[ for n > 1. 
However, the knowledge of g on the interval [0, l/2[ is not sufficient to determine g on the 
next interval [1/2, 2/3 [, as the derivative of g is not continuous at a; = 1/2. For instance, for 
/? > 0, one has g' (^ ) = 0, whereas it can be shown that g'{l/2 + e) ~ e~'^\og{e). As the 
differential equation in Eq. ( pi]) is of the second order type, one of the integration constants 
remains unknown, the other being determined by the continuity condition at x = 1/2. We 
are thus left with a classical shooting problem, where g'{2/3) will be fixed by asking that 
the distribution vanishes at x = 1. 

Let us make this point clearer in the case /? = 0, for which the first iteration can be 
explicitly performed starting from g{x) = a = g{0) for x G [0, l/2[, where the constant a 
is given in Eq. (18). This constant can be eventually calculated once the full distribution 
is known up to this overall constant, as it will ensure the proper normalization of the 
distribution /(x) = x~^g{x). For x G [1/2,2/3], we get. 



g{Q) Ax 



^ + ^-2-2 



X 




X 



2(2x - 1) 



+ c 



(23) 



where c = — ^- ^ ^,L ' . Note that in this case, we indeed find that g'{l/2+e) ~ log(£), leading 
to an infinite derivative at x = (1/2)+. This form for (?, valid on the interval [1/2, 2/3], the 
result of Eq. (pT]) and the continuity condition for g and g' at u„ (for n> 1) leads to the full 
determination of g. Then, a proper choice of the constant c ensures that /(I) = 0. 

In fig. 1, g{x) (equal to /(x) for /5 = 0) is shown for different values of c. For the 
optimal choice for c, it coincides perfectly with the distribution obtained by directly iterating 
Eq.(|3D. 

In fig. 2, g[x) = x^/(x) is shown for /3 = —1/2 (weak disorder) and /? = 1/2 (strong 
disorder). The small x behavior for / is confirmed, and the distributions are again in perfect 
agreement with the numerical integration of Eq. (|2l|) . 



C. Lifshitz tail at the edge of the spectrum 

Note that all these distributions seem to vanish well before Xmax = 1- We will show below 
that this is not the case and that this apparent behavior can be accounted by the fact that. 



/(I 



-i/e 



when e -^ 0. 



(24) 



We shall see in section III.E that the apparent multifractal properties of the Xj's observed 
in [^ are partly due to this phenomenon. 

Let us give a qualitative justification of Eq. (pi]). Taking x = 1 — e, and expecting a 
very fast decay oi g aX x = 1, the two most singular terms in Eq. (^TJ) should be g" and the 
RHS. It can then be checked that the ansatz g{l — e) ^ e^°''^ is solution of Eq. (|21| ) (up to 
subleading multiplicative logarithmic terms) if one takes a = 4. 

Still, actually solving Eq. (^), even in the limit x ^ 1, remains a formidable task and 
the previous argument should be taken with care. However, we can justify rigorously the 
fast decaying tail of / at x = 1, finding a result fully compatible with Eq. (pi]). 

Consider P{e) = Sl„^f{x) dx, the probability to have Xi > 1 — e. In a way similar to 
that leading to Eq. (^OD , we find. 



.1 .1/2 .1/2 

P{e) = dy ^^^dq dpf{y)p{p)p{q), 25 

Jl-he) J l-he)) 2v J 1-e) 1-pv) 



where h{e) = £:/(! — e). For small e, we expect P{6) to be very small as p and q are to be 
taken close to 1/2 and, simultaneously, y must be close to 1 (see Eq. (|13D). In the vicinity of 
(p, q) = 1/2, the distribution p is smooth and roughly constant (see Eq. (|^)). We thus get: 

P{e) ~ 4(1 - Pf f' dy f'^' dq C'^ dpfiy), (26) 



ll-h{e) J{l^h{s))/2y J{l-e){l-qy) 

ll-h(e) 



f ^dyf{y)[l-y-h{e)f. (27) 

Jl—Me) 



From now, we use the symbol ~ in its true mathematical sense, such that a{e) ~ h{£) means 
that a{e)/h{£) —>■ 1, when e —>■ 0. For sufficiently small e, this leads to the following upper 
bound: 

P{e)<a^l' dyf{y)[l-y-h{e)]\ (28) 

Jl-h{e) 

valid for any constant a+ > *■ ~ -* . We then get, 

»l-e „ /-l 

-h{e) 

where we have used the fact that, 



P{e) < a+ /'" dyfiy) [s - h{e)f + a+ /' dyf{y)h{sf, (29) 



[l-y-hie)f<[E-hie)fr^e\ for y e [1 - h{e),l - e], (30) 

[l-y-h{e)f<hieYr^e', for ye[l-e,l]. (31) 

Thus, there exists two constants ai and 02, such that, 

P{e) < aie^P{h{e)) + a2e^P{e). (32) 

Finally, this last inequality shows that for e sufficiently small, there exists a constant A^ > 
such that, 

P{e) < A+6^P{h{e)). (33) 



On the other hand, using again Eq. (pT]) and choosing a sufficiently small S to be determined 
below, we have. 

Pie) > ^ lUie)ii+s) dyf{y)[l - y - hie)f , (34) 

>A_P{h{e){l + 6))e^6^ (35) 

again valid for any constant A^ < ^ ~^' , for small enough e. In the following, we take 
5 = 6°", and will fix the constraint on a later. 

Let us now start from a small enough e such that the inequalities in Eq. (|33| ) and Eq. ( p5D 
hold. Taking 6^=6^= e, we then define e^ and £" by the recursion relations. 

Both sequences go to £ = 0. These recursion relations can also be rewritten, 

-^-^ = 1, (37) 

^-^ = l + 0(e'^J-'^). (38) 

Thus, for any a > 1, both sequences are equivalent to n^^ (by applying Cesaro's mean 
theorem): 

el^e-^ ^. (39) 

Then, by iterating the recursion relations of Eq. (^,^), we get: 

n-l 

log(P(4)) < 4 ^ log(4) + nlog(A+) + log(P(e)) ~ -4nlog(n), (40) 

fc=0 
n-l 

log(P(£;)) > 2(1 + a) ^ log(£;) + nlog(A_) + log(P(£)) ~ -2(1 + a)nlog(n). (41) 

fc=0 

Finally, using Eq. ( |37| , |38D and the fact that P{e) is a continuous and increasing function of 
£, and as a can be arbitrary close to 1, we get that for any arbitrary small r/ > 0, there 
exists £ > 0, such that for any < £ < e, 

_ 111^ iog(5) < - log(P(e)) < _lil±:Zl log(e), (42) 

which leads to our final result: 

-log(P(£))~-log(/(l-£))~-^log(£), when e ^ 0, (43) 

which is a more precise and rigorous statement than that of Eq. (0). It can also be shown 
that the subleading terms in Eq. (^) are a priori oi order log(log(l/£:))/£:. In practice, these 
strong subleading corrections and the very fast decay of the distribution near a; = 1 make 
the quantitative numerical confirmation of Eq. (^31) quite difficult. 

8 



D. Escape probability distribution: symmetric case 

In this section we are interested in the symmetric version of our model for which Qi = Pi-i. 
The approach is shghtly different from that of the previous case but is definitively in the same 
spirit. As a consequence, less attention will be attached to rigorous arguments, although 
they can be adapted without any difficulty to this problem. 

The map now becomes, 

x^ = z ^ , (44) 

l-pi_lXj„i 

which shows that the novel variable yi = PiXi satisfies the recursion, 

m = T-^- (45) 

1 - Vi-i 

u = p^ has a distribution cr(n) satisfying, 

.1/2 



Ci[U] 



dw piw)5iu - w^) = 2-/^(1 - P)u~^'+''^/^9{u)e{l/4: - u). (46) 



Using Eq. (^5]), we find that y is always in [0, 1/2], and that for y G [0, 1/4], the probability 
distribution F{y) of y satisfies, 

F{y)=2~^{l-(3)y~^'^^y' f'^' F{y'){l - y'f-^^' dy' , (47) 

Jo 

that is a pure power law behavior. For y > 1/4, and proceeding along the same line as in the 
non symmetric case, we find that F{y) satisfies the following non local differential equation: 



2y ^^^ 8y3 I Ay^ 



F'(y) + ^F(y) = ^F I 1 - - 1 . (48) 



Again, the knowledge of F{y) on the interval [0,1/4], permits the determination of the 
distribution on the next interval [1/4, 1/3], and by recursion on each of the intervals of the 
form [un,Un+i], with Un = 2{^- 

Now, let us analyze the behavior of P{s) = JiL-e F{y) dy, for small e. Defining h{e) = 
s/{l — 2e), we easily find that there exists two constants c and C (which can be actually 
determined) such that, 

P{e) r^c I ' {h{e) - z)F{l/2 -z)dzr^ Ce^P{h{e)). (49) 

JO 

The last estimate is obtained using the same type of inequalities as in the non symmetric 
case. Again defining, Sn+i = £n/(l + 2£:„), and using the fact that e„ ~ {2n)~^, we find that, 

log(P(en)) ~ 2 Y: log(e.) ~ ^^^, (50) 

fc=l ^" 

which finally shows that. 



log(P(e)) ~ log(F(l/2 - e)) ~ ^. (51) 

From the knowledge of the properties of F{y), we can access to that of f{x), the stationary 
distribution of the Xj's, using the relation, 

fix) = fj^ dpfj^ dyp{p)F{y)6 (x - ^] . (52) 

Let us exhibit the main properties of / which will be completely similar as what was obtained 
in the non symmetric version of the model. For x G [0, 1/2], Eq. (^) leads to, 

fix) = 2^-^(1 - P)x-^ f'^' F(y)(l - y)'-^ dy, (53) 

Jo 

which shows that f{x) is simply proportional to x~^ as in the non symmetric case. For 
X e [1/2, 1], we obtain, 

g'ix) = ix^fix))' = ^F (l - ^) . (54) 

Using Eq. (pl|), we finally conclude that, 

log(/(l -e)) ~ log (f (l - ,rr^]] ~ -log(£), (55) 



that is a similar behavior as was found in the non symmetric case (see Eq. (^31) ) up to the 
factor 4 which is replaced by a factor 2. The physical interpretation of this smaller coefficient 
is quite clear: in the non symmetric case, for Xi to be close to 1, one must have Xi-i, pi and 
Qi close to their maximal value. In the symmetric case, for Xi to be close to 1, only Xj„i and 
Pi must be close to their maximal value, as g, = pj_i is automatically close to 1/2 as Xj_i is 
close to 1. This extra constraint explains why f{x) decays faster in the non symmetric case 
than in the symmetric case which is confirmed numerically on fig. 3. 

Note that the non symmetric case could have been treated by the same technics as in 
this subsection by replacing the distribution a oipj by the distribution a oi piQi+i. 

Finally, we can conclude that up to a few irrelevant details, the symmetric and non 
symmetric models seem to share exactly the same properties. This is apparently surprising 
as the corresponding models without trapping are drastically different 0. This intriguing 
property is confirmed and explained physically in section IV.D. 

E. Explanation of the apparent multifractality 

In the preceding subsections, we have obtained a puzzling numerical result (see fig. 1-3): 
although we have shown that the maximum value Xmax = 1 must be attained, the numerical 
maximum value effectively obtained after 2 x 10^ iterations of the map Eq. (13) is typically 
of order Xmax.eff ~ 0.73 ~ 78, for the three values of (3 actually tested in the non symmetric 
case. This apparent paradox can be explained by the sharp decay of the distribution f{x). 
Xmax.cff can be estimated by considering that after A^ iterations of the map, 

10 



f{x)dxr^N-\ (56) 

which is the smallest non zero value that this integral can take (the integrated distribution 
increasing by elementary steps of height A^~^). Using Eq. (0), we find that e = 1 — Xmax,eff 
must satisfy, 

-!2i^~ilog(iV). (57) 

If we now take N = 2 x 10^, the above estimate gives Xmax,eflF ~ 0.75, in fair agreement with 
the observed range of numerical effective values for Xmax- 

In the limit of very large N, Eq. (^) also leads to the leading order estimate, 

,. 41og(log(iV)) 
l-a:^ax,efr log(iv) ' ^^^) 

which goes to zero very slowly. 

In 0, the fact that Xmax = 1 was not recognized. The authors actually computed the 
multifractal distribution of the Xj's in the interval [O,a;max,off]- This interval was cut into n 
equal length intervals, and pj was defined as the fraction of the total number of the Xj's 
which belongs to the j^^ interval. One then defines [|1^, 



n 



Z{q,n)=Y.plr.n-^^^\ (59) 



i=l 



the last equivalent defining the scaling exponent r(g) associated to the q^^ moment. Let us 
first derive the exact expression of r(g) for a given choice of /5 > (strong disorder), and 
taking the actual value of Xmax = 1- For < g < 1//5, the function f\xY is integrable, so 
that, 

Z{q,n) = n-(''-i) x -Y^npi)" ~ n^^'?"^) / fixfdx. (60) 

n .^^ 7o 

This shows that in this regime, r(g) = g — 1. For q > 1/(3, the function /(x)"^ is no longer 
integrable due to the power-law divergence at x = 0: 

Z{q, n) = n-^"-^) x - J^i^PiY ~ ^"^'"'^ [' fi^Y dx ~ n^^^^^)". (61) 

In this regime, we thus find r(g) = (1 — j3)q. Morover, for q = 1//3, we find, 

Z{q = 1/(3, n) ~ n-(^/^-i) /' /(x)^/^ dx ~ n-^^/^"^) log(n), (62) 

Jl/n 

which shows that in the region q ^ 1/(3, the numerical determination of r{q) will be strongly 
affected by a logarithmic slow cross-over. Finally, strictly speaking, r(q) is not defined for 
negative q, due to the essential singularity at x^ax = 1- But if one computes the multifractal 
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scaling exponents by restraining the study on the interval [0,Xmax,eflf] (as in 0), we then 
recover r(g) = g — 1, for negative values of q as well. 

The two linear regimes for r(g) are clearly visible on the fig. 4 of |p, with the predicted 
slopes and transition point q = 1/(3. Of course, numerically, the change in slope (from 
T'(g) = 1 to r'(g) = 1 — /3) is found to be smooth, partly due to the logarithmic correction 
around q = 1/(3. Then, the spectrum of singularities defined as the Legendre transform of 
r{q) 0, 

s(a) = ag — r(g), a = r'(g) (63) 

apparently yields a non trivial spectrum, whereas the actual one is concentrated on two 
points: a = 1, with a support of fractal dimension s(l) = 1, reflecting that for almost all 
values of x the distribution /(x) is actually continuous, and a = 1 — (3, with a support 
of fractal dimension s(l — /5) = 0, which simply results from the x = singularity of the 
distribution /(x). 

The moral that we can draw from this is that one must be very careful when dealing with 
multifractal analysis, especially if there is no special reason to expect a multifractal spec- 
trum. Similar problems were encountered by the author of |]T^, who obtained an apparent 



multifractal spectrum in a model which can actually be solved exactly |T^, and for which 
it can be shown that the true multifractal spectrum is of the same type as above. Similar 
doubts can be raised on the findings of multifractal spectra in certain biological systems or 



in the field of finance [T^, where simple power-law distributions can lead to such apparent 
behaviors. 



IV. LARGE TIME BEHAVIOR OF THE SURVIVAL RETURN PROBABILITY 

A. General results 



Following the tracks of [|T^, let us evaluate the survival return, or more exactly, its 



discrete time Laplace transform. We consider the symmetric model but the non symmetric 
case can be treated in the same spirit, leading to exactly the same results. 



Thus consider, 



Po.H^J 



g^cMl^L. (64) 



It can be shown [|T0| that pQ^i{uj) satisfies the equation (see Eq. (|^), 

Po,i(^) = ti'i-i,jPo,j-i(^) + Wi^i+iPo^i+i{u) + (1 + cj)"^(5i,o- (65) 

Then, the variables (t)t{!^) and (j)~{uj), defined respectively for i > and z < by, 

(/),^(u;) = ■ -— , (/),• Ujj) = ■ -— , (66) 

satisfy the same recursion, reminiscent of that of Eq. (|l^). For instance, 
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/^?-l 



l-0*+l(^) 1 + ^ 

with a similar equation for i < 0. It can then be shown that, 



(Po,o(^)) 



n„ 



)n„ 



<>+ - 



(67) 



(68) 



where 11^(0), is the expected stationary distribution of 0+ and 0". 

From now, we follow the method of section III to evaluate the behavior of the distribution 
n^(0) close to = 0max, and for small uo (as we are mainly interested in the large time 
behavior of Ps(t)). From Eq. (|67|) , 0max(^) can be easily calculated: 



(69) 



</>ma.(^) = ^ (l - Vl-(1+^)-') = ^ - /I + ^(^'^') 



Moreover, using again Eq. (|67|) , we find that 11(^(0), satisfies the following self-consistent 
equation: 



n.(0) 



1/4 /•</>max(t<j) / 

(ix / (i0'o"(x)nu;(0')5 

JO \ 



X 



:i+'^)^(i 



(70) 



where the distribution a{x) oi x = yU^(l + uY is given by Eq. (^). 

Eq. (|67|) can be first solved for < ^.-^ );^.2 , leading to a pure power-law behavior (as 
encountered in section III): 



n,(0) = 2-^(1 -f3){i + cu)-(^+^)0"^'+''^/' 



<^max{tij) 



n.(0) 



0/)(l+/3)/2- 



For > ^^^^.2 ; one can differentiate Eq. (|70|) (noticing that 0' G 
leads to the following differential equation for 11^^(0): 



1 



4(l+a;)2 



0(^+^)/^n.(0) 



1-/3 



(l + ^)-(5+/3)/2^-{3-/3)/2nJl 



1 



4(1 + cj) 



(71) 
which 

(72) 



The leading asymptotics of Iluj{(p) close to = 0max(^) can be calculated rigorously adapting 
the method of section III.C. Here, we first derive this result by a less rigorous method 
already mentioned in section III. A, consisting in keeping only the most singular terms in 
the differential equation Eq. ([7^), leading to: 



""(^'-^""(i-iirT;;)^) 



' ^n^ 0max(^) - 



1 - 0max(a;) 
0max(t^) 



(73) 



where A is a computable positive constant, and the explicit equation for 0max('^) was used. 
One can try an ansatz form for na;(0) which satisfies this equation up to multiplicative 
logarithmic terms. We find: 



iog[n. 



-uj V Y'lnax 



.00] 



e)] 



\og\e) 
2 log(r) 



with 



- 0max(^) 
0max(c^) 



(74) 
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B. Tail of the distribution 

In fact, this complicated ansatz was originally found by applying a similar method as 
that of section III.C, that we present now. 

Again, let us define Pu){,£) = ^^^^^^^_^Il^{(f)) d(f), which satisfies: 

/■1/4 /-^maxCl^) ( X \ 

P^{e) = j^ dx j^ d<P a{x)Yi^{4>)e f +e- <pra.^{uj) \ . (75) 

Once we express the actual domain of integration by imposing that the argument of the 
6 function be positive, we find that the variable x remains very close to 1/4 where the 
distribution (t{x) is essentially constant. Exploiting this fact, we find after a few elementary 
manipulations that, 

pre 

P^{e)^C P^{u)du, (76) 

Jo 

with r = "'^'"^ \ > 1) ^-iid C = ."}(■. ■ This equation is a rigorous integrated version of 
Eq. (ffSl). Note that for small uo, we have, 

r = 1 + 2^2^ + 0{uj'^'^). (77) 

We can now proceed in the same spirit as we did in section III.C, and find exact inequalities 

for Pt^{e). For any c+ > C, and sufficiently small e\ 

Pu.{e) < c+ /r PM du + c+ /o^ PM du, (78) 

< c+ir - l)EP^ire) + c+eP^{e), (79) 

which finally leads, for small enough e, to the existence of a constant C+ of order 0(1), such 
that, 

Pu.{e)<C+{r-l)eP^{re). (80) 

On the other hand, for any r' < 1 close to 1, to be specified later, we can write that, 

P^{e)>C-r,f^^P^{u)du, (81) 

>C-{l-r')eP^{r're), (82) 



Let us now start from a small enough 6^ = 6 such that the inequalities of Eq. (^,^) hold, 
and define, 

4 = ^"'4-1 = r-'^e, (83) 

e- = {r'r)-^e-_, = {r'rr^e. (84) 

In the following, we will choose r' such that r'r > 1, so that e~ -^ 0, when n — >■ +oo. 
By iterating the recursion inequalities we get. 
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n-1 



log(P^(4)) > -nlog(C+(r - 1)) - E log(4) - log(il.(£)), (85) 



fc=0 
,2 



?2 

~ -nlog(C+(r - 1)) + y log(r), (86) 

where, after using n ~ — "^ /^-, , the last hne can also be written, 

loS^(4) , Mg+(r-l)) . +. . . 

Similarly, 

- log(P.(£;)) < -nlog(C_(l - r')) - E^=o' log(efc ) - log(P.(e)), (88) 

2 

where the last line can also be written under the form 

log2(£;) , log(C_(l-r')) 



nlog(C_(l-r')) + ^log(r'r), (89) 



21og(r'r)^ log(r'r) ^°S(^n)- (90) 

As r' can be arbitrarily close to 1, we thus find that the leading order of Eq. ( |7iD is exactly 
recovered. 

Now, we restrict ourselves to the case of small uj, and analyse the effect of the subleading 
term. To be specific, let us take 1 — r' = cu^''^^^, with 5 arbitrarily small, such that the 
condition r'r > 1 remains satisfied (see Eq. ([77D). Our final results is that P^{e) < Puj{e) < 
P+{e), with, 

-log(P.^(.))^^ + ^log(.), (91) 



- log(P„-(.))~ i^ + (1 + 2^)M^log(.), (92) 



which strongly suggests that. 



log(n^(0^ax(cu) - e)) ~ - log(P.(e)) ~ l"g(^) l^i^^^^ (93) 



Note that the case u = falls exactly in the class of problem studied in section III.C 
leading to the exact asymptotics (with 0max(O) = 1/2): 

log(n^=o(l/2 - e)) ~ log(P^=o(^)) ~ ^^^. (94) 
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C. Survival return probability 

The results of this section will not rely on as firm mathematical grounds as that of the 
preceding sections, but will appear to be quite reasonable. 



The comparison of Eq. (^) and Eq. (Q) suggests that the asymptotics of Eq. (^) should 
be correct at least up to e of order ^/cJ ~ 1/2 — ^maxlto") (or more exactly nw^)- Using 
Eq. (|68|), we find that there should be a singular contribution in (po,o(^)) of order (see 
Eq.(El), 



[(Po,o(0)) - (Po,o(^))].,,, - r ^ ^ n_o(a:) dx ~ exp ("Hli^] , (95) 



J Sing 



where the last estimate comes from Eq. (|9J). In principle, the leading singular correction to 
Eq. ( p5D should come from the contribution of 11^ to the integral of Eq. (p8|) on the interval 
[0, (/'max] and should be of order exp f — ^^ °^ )■ 

A contribution of the form exp(— c| log(ti;)|" 1^°") in a Laplace transform generally origi- 
nates from a large time decay of the form (as found by a steepest descent type of argument). 



exp ( -CtTT^ log w^ (t) J . (96) 

If we take the logarithmic corrections in Eq. (|95D seriously, we thus find: 

(po,o(t))~exp(-Cti/3log2/=^(t)), (97) 

in disagreement with the heuristic argument given in [^ (see Eq. (^). 

D. Exact bound for the survival return probability 

In this subsection, we give an exact lower bound for the number of surviving particles 
using an argument which can be easily generalized to obtain a similar bound for (po,o(^))- 
The resulting bound is fully consistent with Eq. (|97D and contradicts the heuristic estimate 
of[g(seeEq.(|)). 

Consider the symmetric model (wj^+i = Wj+i,i). The probability to have a region of L 
sites on which all ^1^2, • • • , 'W^l-2,l-i > 1/2 — ^ is. 



VAe) 



L-2 



l-{l-2eY-^ ~ [2(1 - /3)£]^. (98) 



If we were interested in the non symmetric model, the corresponding probability to have 

wi,2, ^2,1 . . . , wl-2,l-i, wl-i,l-2 > 1/2- e would be simply. 



2(L-2) 



VL{e) = [1 - (1 - 2ey-^\ ~ [2(1 - /?)£]^^ (99) 

and the rest of the argument would be essentially identical. On such a region, the density 
cannot decay faster than that of the following problem where we consider Wi^2, ■ ■ ■ , wl-i,l = 
1/2 — e, with fully absorbing boundary conditions (u^o,i = Wi^ = w^^i^l = "U^l.l-i = 0). 
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This simple property can be shown by induction using the master equation for Pi{t), the 
density at site i: 



Pi{t + 1) = Wj+i,iPi+i(t) + Wi_i^iPi_i{t), pi{t = 0) = 1. 



(100) 



This simpler problem can be solved exactly (on the lattice or in the continuum), leading to 
the following bound for the average decay of the total density: 



riiit) > ^-exp 



TT' 



2lJ 



2et 



;ioi) 



Finally, we find that the total density n{t) is bounded for any L and e by the following exact 
lower bound: 



n{t) > Cocxp 



2L2 



2et + Llog{2{l-f3)e) 



(102) 



where Cq is a positive constant. We can now take the maximum of this lower bound over 
L and e. Expressing this condition, we get the following optimal values for L and e defined 
imphcitely by. 



L = 2et 



-kH 



log(2(l-/3)£) 



1/3 



21og(t) 



nl/3 



(103) 



where the last estimate is valid for large time. 
Finally, we have shown that, 

n(t) > exp-5(t), with S{t) 



SttH 



1/3 



log^/^(t), 



(104) 



in perfect agreement with the above analytical argument. The physical interpretation of this 
result is that surviving particles are living in large regions where the Wi^i±i are very close to 
1/2, and annihilates with a large probability outside these regions. This explains why the 
symmetry of the hopping probabilities is irrelevant and why the results is essentially similar 
as that of perfectly diffusing particles with randomly distributed perfect traps. The log^'^(t) 
corrections is due to the fact that for large time, the Wi^i±i in regions where the surviving 
particles stand must be closer and closer to 1/2, with an allowed fluctuation decreasing as 
€ ~ t~^/^. Moreover, within these large regions there is an extra probability leaking per site 
of order e. Note finally, that if < 7 < 1, the argument can be repeated with et replaced 
by (1 — 7)£:t in Eq. ( |101| ), leading to the same decaying behavior. 

Let us now exhibit the fiaw(s) in the argument given in [^, leading to Eq. (|^). On a 



large region of size L with Wi^i±i > 1/2 — e, the authors of [|T0| estimated the probability 
decay as. 



niit) ~ exp 



LV^ 



eL^/^ 



(105) 



where u is the effective diffusion exponent defined in Eq. (^) and below. This estimate is 
to be compared to our Eq. (jlOll) . The first term is supposed to represent the probability 
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decay due to the absorption of particles at the boundaries of the considered region. It is 
not correct to consider that particles in this region display anomalous diffusion as the Wi^i±i 
are in fact very close to 1/2. This fact is confirmed by the exact bound Eq. ( |101| ). The 
second term in the exponential is supposed to reflect the fact that there is a small trapping 



probability (of order e) on each site of the considered region. The authors of |T0[ assumed 
that time can be replaced by L^^'^. This is obviously wrong as the decay due to this small 
trapping probability explicitly depends on the time spent in the region but not on its size. 
Finally, the authors did not realize that in the final expression that they obtained, e (and 
not only L) should also be treated as a variational parameter. 

Note that for intermediate times, we expect that the density should decay as, 

n{t) ~ exp(-cA^(t)) ~ exp{-Ct^/^), (106) 

where N{t) is the number of different sites visited by a random walker. This phenomenon 
is also well-known for the Donsker-Varadhan problem, for which this behavior can actually 
dominate the numerically accessible time regime |0,^. 

Finally, the generalization to higher dimensions of this model is straightforward. On a 
periodic lattice of coordination number z, the hopping probabilities are bounded by z~^ and 
particles disappear if they do not move. The above argument suggests that, 

n{t) ~ exp (-Ct3T2 logdT2 (t)) , (107) 



V. CONCLUSION 

In this paper, we have considered a model where the trapping probabilities are strongly 
correlated with the hopping probabilities of the walker. We have shown that the escape 
probabilities have a well defined distribution which has been analyzed in great detail in 
section III. We have also explained why this quantity display an apparent multifractal 
distribution. In section IV, we generalized our approach to the study of the survival return 
probability distribution. We deduced from this exact analysis that this survival probability 
(and the density) decays as exp(— Ct^/^log^'^(t)). To support this result, we have obtained 
an exact bound for the density which even reproduces the correct power of the logarithmic 
correction. Moreover, we have explained the independence of the model properties with 
respect to the disorder strength /?, the trapping rate 7 > and, more surprisingly, the 
symmetry of the hopping probabilities. 

A challenging problem is the understanding of the diffusion properties in this model. 
The fact that the effective spreading of the survivors is apparently faster than diffusive (see 
Eq. (ID) remains to be explained. 
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FIG. 1. In the case /? = (/(x) = g{x)), we plot the distribution obtained by iterating Eq. (13) 
2 X 10^ times (N.S.), and the solution of Eq. (^ij) obtained by imposing /(I) = 0, which leads to 
c = 0.0678. . . (see text). The two curves are indiscernable. We also plot the normalized solution 
of Eq. (0) for c = 0.04 and c = 0.10. 
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FIG. 2. We plot g{x) = x^f{x) for /3 = —1/2,0, 1/2. In each case, the distribution is obtained 
by iterating 2 x 10^ times the map of Eq. (p!3|). The agreement with the numerical solution of 



Eq. (21) is perfect. 
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FIG. 3. We plot g{x) = x^ f{x), for /3 = 1/2, and after 2 x 10^ iterations of the symmetric and 
non symmetric maps. As explained in the text, the distribution decays faster in the non symmetric 
case. 
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